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ABSTRACT 

We explore the morphology of Type la supernova remnants (SNRs) using three- 
dimensional hydrodynamics modeling and an exponential density profile. Our model 
distinguishes ejecta from the interstellar medium (ISM), and tracks the ionization age 
of shocked ejecta, both of which allow for additional analysis of the simulated rem- 
nants. We also include the adiabatic index 7 as a free parameter, which affects the 
compressibility of the fluid and emulates the efliciency of cosmic ray acceleration by 
shock fronts. In addition to generating 3-D images of the simulations, we compute 
line-of-sight projections through the remnants for comparison against observations of 
Tycho's SNR and SN 1006. We find that several features observed in these two rem- 
nants, such as the separation between the fluid discontinuities and the presence of 
ejecta knots ahead of the forward shock, can be generated by smooth ejecta without 
any initial dumpiness. Our results are consistent with SN 1006 being dynamically 
younger than Tycho's SNR, and more eSiciently accelerating cosmic rays at its for- 
ward shock. We conclude that dumpiness is not a necessary condition to reproduce 
many observed features of Type la supernova remnants, particularly the radial proflles 
and the fleecy emission from ejecta at the central region of both remnants. 

Key words: Hydrodynamics - instabilities - ISM: supernova remnants. 



1 INTRODUCTION 

SN 1572, visible to the unaided eye on Earth, is now linked 
to the name of its most accurate observer, Tycho Brahe. 
Its remnant was rediscovered in radio frequencies almost 
four centuries later (Hanbury-Brown & Hazard 1952; Bald- 
win & Edge 1957), and has since been observed throughout 
the electromagnetic spectrum. SN 1006 was probably signif- 
icantly brighter than SN 1572 to observers (Krause et al. 
2008; Winkler, Gupta & Long 2003) around the globe when 
it first appeared, but its remnant was not identified until 
Gardner & Milne (1965). Both supernova remnants (here- 
after SNRs) are resolvable to several hundred pixels across 
their diameters with the Chandra X-ray Observatory, and 
confirmation that the original supernovae were of type la 
(Baade 1945; Krause et al. 2008; Wu et al. 1983) make them 
ideal laboratories to test theories of supernova and remnant 
evolution. 

Detailed observations have shown both Tycho's SNR 
and SN 1006 to be rich in structure. Of the numerous 
features observed in Tycho (Seward, Gorenstein & Tucker 
1983; Dickel, van Breugel & Strom 1991; Reynoso et al. 
1997; Warren et al. 2005) and SN 1006 (see Cassam-Chenai' 
et al. (2008) and references therein), three aspects of the 
remnants' morphology are relevant to this paper. One is 



the nature of the clumpy, "fleecy" thermal emission located 
throughout Tycho and SN 1006, long since identified as 
ejecta in Tycho's case (Hamilton, Sarazin & Szymkowiak 
1986) and that could be caused by the Rayleigh- Taylor in- 
stability. Another is the limb brightening effect, which can 
be used to identify the reverse shock but which is less pro- 
nounced than expected, or even absent, in many places 
around Tycho and SN 1006. The third is the radial structure 
of the remnants; that is, identifying the radial locations of 
the forward shock, contact discontinuity, and reverse shock 
and their projected positions. 

It has long been argued that SNRs would be subject 
to fluid instabilities such as the Rayleigh- Taylor (Gull 1973) 
and the Kelvin-Helmholtz shear instability (EYyxell et al. 
1991; Chevalier et al. 1992). Swept-up interstellar mate- 
rial (ISM) decelerates the denser ejecta behind the con- 
tact discontinuity, and fluctuations within the ejecta or ISM 
seed the Rayleigh- Taylor instability. Early treatment of this 
problem was either analytical or limited to one-dimensional 
numerical simulations looking for Rayleigh- Taylor-unstable 
zones in a radial proflle (Gull 1973; Dickel et al. 1989). Non- 
linearity and asymmetry limit the effectiveness of either ap- 
proach, but advances in computing have since allowed SNRs 
and relevant instabilities to be modeled in multiple dimen- 
sions (Hachisu et al. 1990; Pryxell et al. 1991; Chevalier et 
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al. 1992). These results were expanded upon by other au- 
thors for a variety of other scenarios, such as ejecta clump- 
ing (Orlando et al. 2012), ovcrdcnsc cjccta "bullets" (Wang 
& Chevalier 2001), underdense ejecta "bubbles" (Blondin, 
Borkowski & Reynolds 2001), efficient particle acceleration 
at the shock fronts (Blondin & Ellison 2001; Ferrand et al. 
2010), SNR interactions with a moving ISM (Velazquez et 
al. 2006), or asymmetries in the progenitor star and its en- 
vironment (Blondin, Lundqvist & Chevalier 1996; Vigh et 
al. 2011). 

SNRs expand and interact in three dimensions, but are 
observed in only two. Thus any comparison against obser- 
vations must include some projection of the data, leading to 
projection effects such as limb brightening or overemphasis 
of irregular surfaces. Ferrand et al. (2010) present a three- 
dimensional simulation of one octant of a remnant, including 
line-of-sight projections of radiating material. Images both 
with and without emission from shocked ISM can be found 
in that paper, and show developed ejecta structures similar 
to those in Tycho. They do not, however, accurately capture 
the weak limb brightening effect observed by Warren et al. 
(2005); additionally, they find the reverse shock closer to the 
forward shock than observations suggest is the case for Ty- 
cho. The work of Orlando et al. (2012), performed over 471- 
steradians, uses volume renderings of the evolved remnants 
to match the radial structure of SN 1006. Inclusion of ejecta 
clumps alters the corrugation of the reverse shock, with the 
amount of disruption increasing as the size of the clumps in- 
creases. It also allows Rayleigh- Taylor structures generated 
by these clumps to reach and perturb the forward shock. 

The primary question addressed by this paper is 
whether fluid instabilities alone are sufficient to explain the 
morphology of Tycho's SNR and SN 1006 as observed by 
Chandra. All three features of interest in these SNRs the 
proximity of the contact discontinuity to the forward shock 
in projection, the fleecy nature of the central emission, and 
the indistinct limb brightening in some locations - would 
be affected by fluid instabilities at the contact discontinuity. 
Projection of a three-dimensional shell structure to a two- 
dimensional surface overemphasizes deviations from spheri- 
cal symmetry, such as Rayleigh- Taylor fingers; this leads to 
errors in estimating the locations of the contact disconti- 
nuity and reverse shock. In particular, it has recently been 
proposed (Orlando et al. 2012) that clumpy ejecta is nec- 
essary to explain these features. We present in this paper 
complete three-dimensional simulations of a type la super- 
nova remnant, covering the full angular range to allow unre- 
stricted growth of Rayleigh-Taylor instabilities in all direc- 
tions, without preferential treatment of certain wavenum- 
bers. 

In §2 of this article, we explain the model used in the 
simulations. In §3 we describe the initial conditions of the 
simulations and present results from the runs. In §4 the 
three-dimensional simulations are projected into a plane and 
compared against observations of Tycho's SNR and SN 1006. 
We discuss features of the projections in §5, in particular 
our estimation of the dynamical age of the two remnants. 
We conclude in §6. 



2 THE EJECTA PROFILE 

Dwarkadas & Chevalier (1998) compared several models 
for the ejecta profile of a type la supernova - constant, 
power law, and a new exponential profile for ejecta density 
- against several hydrodynamical models for type la explo- 
sion mechanisms and early spectral evolution. They found 
that the exponential density profile fit the models under 
consideration to a higher degree of accuracy than did the 
constant or power law profiles. In particular, the authors 
noted a good fit with the successful W7 model of Nomoto, 
Thielemann & Yokoi (1984) - see Figure 1 of Dwarkadas & 
ChevaUer (1998) - and found that the profile began steep 
but flattened out over time. We adopt the exponential model 
developed in that paper, which uses a spherically symmetric 
density profile 



PSN = At e 



(1) 



Assuming that r = vt aX t = to, the constants A and Ve are 
found to be 
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where Mch 1.4 and £51 are, respectively, the Chan- 
drasekhar mass and the supernova energy in units of 

lO''^ ergs. 

We use the same scaling as Dwarkadas & Chevalier 
(1998): 
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with no = pam/(2.34 x 10^^* g) being the number density of 
the interstellar medium, assuming a 10:1 H:He ratio. Pres- 
sure values were scaled to the ram pressure PamV'^. This 
scaling is not unique; McKee & Truelove (1995) adopt a dif- 
ferent scaling in their treatment of SNRs. For this paper 
we adopt the convention that primed lowercase letters de- 
note the real, physical quantities, while unprimed lowercase 
letters are used for their scaled counterparts, e.g. t = t' /T' . 

Though Dwarkadas & Chevalier (1998) considered mul- 
tiple different ejecta profiles in their one-dimensional study, 
we use only their exponential profile here. Other work (Fer- 
rand et al. 2010) has been performed in 3-D using the power 
law profile, against which the results presented in this paper 
can be compared. To our knowledge, no multi-dimensional 
simulations have used the constant profile. 
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3 THE NUMERICAL MODEL 

To run the simulations we used the hydrodynamics code VH- 
1, which solves the Euler equations for fluid flow using the 

picccwisc parabolic method with a Lagrangian remap step - 
see Colella & Woodward (1984) for a thorough description of 
the procedure. The base hydrodynamics code was extended 
with a subroutine to calculate ionization age of ejecta ele- 
ments, which could then be used to generate emission maps 
for comparison against observations. 

In order to save computing time, simulations were initi- 
ated in one dimension, then continued in 3-D at a later time. 
The starting time used in one dimension was one day, or 
t = 2.74 X 10"^ yr/T' = 1.1 x IQ-"^ (Me/Mch)''^''^ El^nl^^ . 
The grid at initialization started at rmin = 5 x 10~^ with 
7"max/''min = 2.4, SO that the forwaxd edge of the grid was 
just beyond where the exponential profile reaches a density 
equal to the ISM. The grid contained 480 zones, for a reso- 
lution of Ar/rmax w 1.2 x 10~^; this was deemed to be an 
acceptable balance between resolution and computing time 
for the three-dimensional runs (a brief discussion of resolu- 
tion effects occurs in the next paragraph, but all 3-D runs 
were performed at this resolution). The initial conditions 
within the ejecta were a radial velocity of v = r/t, density 
set according to the exponential profile, pressure calculated 
as an ideal gas at a temperature of 5 K, and angular ve- 
locities set to 0. The ISM was assumed to be at rest, with 
a density of unity, and cold. The interior radial boundary 
condition was set to the exponential model with no angular 
velocity, and the exterior radial boundary matched the ISM 
described previously. To follow the motion of the shocks over 
many doubling times, the grid tracked the forward shock, 
advancing and expanding by a small amount whenever the 
forward shock moved within six zones of the outer edge of 
the grid. 

Early two-dimensional runs with an effective adiabatic 
index of 7 = 5/3 showed that resolution had a minimal ef- 
fect on the gross structure of the instabilities generated. At 
t = 2.0 in simulations with resolutions between Ar'/rmax ~ 
4.9 X 10^"' and Ar/rmax ~ 2.4 x 10~^, the dominant wave 
mode was consistent across all resolutions, though at higher 
resolutions more fine structure was present. Changing the 
resolution has two effects beyond modifying the power spec- 
trum of Rayleigh- Taylor instabilities. First, increasing the 
resolution would allow for enhanced Kelvin- Helmholtz insta- 
bility formation, which should in turn accelerate the onset of 
the instability saturation period discussed in the next para- 
graph. Second, as the effective adiabatic index decreases, a 
high resolution becomes necessary to resolve the small-scale 
structure associated with highly compressible fluid. 

One constant through all of the two-dimensional runs 
was the presence of three "epochs" regarding instabilities. 
In the first of these, the initial seeds (random perturba- 
tions of either density, radial velocity, or pressure; or a 
long-wavelength perturbation of radial velocity) generated 
Rayleigh- Taylor fingers with a wavelength of a few grid 
zones, which then cascaded into larger structures. The sec- 
ond epoch was a period of instability saturation: multiple 
cycles in which Rayleigh-Taylor fingers appeared, grew to- 
ward the forward shock, experienced shear from the Kelvin- 
Helmholtz instability along their edges, and fell back to- 
wards the reverse shock. An instability "freeze out" divided 



the second and third epochs. After the freeze out, it be- 
came possible to track individual RT structures to the end 
of the simulation, since production of new RT structures 
slowed down and a second cascade to larger structures oc- 
curred. The instability saturation period has been reported 
for both the power law (Chevalier et al. 1992; Kane, Drake & 
Remington 1999; Blondin & Ellison 2001) and exponential 
(Dwarkadas 2000; Wang & Chevalier 2001) models. Because 
of this epoch, we expect that any trace of the initial per- 
turbation, at least for the small magnitudes applied in our 
simulations, is washed away before the fingers freeze out. 

We further investigated the saturation period by chang- 
ing the times at which 1-D runs were expanded to two di- 
mensions and further evolved. One-dimensional runs were 
carried out to t = 3.65 x 10~® (the earliest time at which 
the contact discontinuity appeared as a minimum in the den- 
sity profile rather than as a kink), t = 3.65 x 10^* and 
t = 3.65 X 10"'^ before mapping to two dimensions and 
seeding instabilities. The 2-D run begun at t = 3.65 x 10~* 
achieved the saturation exhibited by the earliest start time. 
The t — 3.65 x 10^^ run never reached the instability sat- 
uration epoch; the first generation of RT fingers were still 
identifiable at the end of the simulation, at t = 2.0. As 
long as the mapping from one to multiple dimensions oc- 
curred early enough, the simulations achieved saturation and 
were largely indistinguishable. We therefore saved computa- 
tion time in the three-dimensional runs by starting them 
at t = 3.65 x 10"'' without any apparent changes to the 
evolution of the remnant. The importance of the instabil- 
ity saturation period to our results bears restating: covering 
four or more decades of time allowed the instabilities to sat- 
urate and wash out traces of the initial perturbations, while 
simulating three decades resulted in clear imprints of the ini- 
tial perturbations in the final state of the ejecta. Orlando et 
al. (2012) used only two decades of expansion time in their 
simulations, from an age of 10 years to an age of 1000 years, 
and their runs with smooth ejecta are still clearly dominated 
by grid effects at late times (note the quadrilateral symme- 
try in figure 6 of that paper, a product of the Cartesian grid 
used in the simulations). 

Given the long-standing evidence in the literature that 
efficient cosmic-ray acceleration at shock fronts can dramat- 
ically effect the eventual morphology of an SNR, three differ- 
ent three-dimensional runs were performed. As an approxi- 
mation for energy loss at shock fronts, we used the procedure 
of Blondin & Ellison (2001), globally adjusting the adiabatic 
index of the fluid to increase its compressibility. The run 
with an adiabatic index 7 = 5/3 (ideal gas with a strong 
shock compression ratio cr = 4) served as a control. To com- 
pare against this we performed simulations with 7 = 4/3 
(<T = 7) and 7 = 6/5 (tr = 11). 

The three-dimensional runs covered the full 4n stera- 
dians, allowing instabilities everywhere in the simulation to 
grow without boundary effects, and preventing grid-induced 
preferred wavenumbers. To avoid geometrical singularities 
at the poles and allow for a more uniform A^ and A0 across 
the grid, we employed specialized Yin- Yang grid of two con- 
gruent parts with angular extent n/2 in 9 and 37r/2 in 4> 
(Kageyama & Sato 2004). A one-dimensional kickstarting 
run was taken out to a scaled time of t = 3.65 x 10"'', then 
swept across the grid to form the basic three-dimensional 
profile, discarding the innermost quarter of the 1-D grid - 
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zones well inside the reverse shock even at the end of the sim- 
ulation. The 3-D runs used a grid resolution of Ar/rmax ~ 
1.2 X 10'^, in keeping with the 1-D kickstart runs. Though 
"cubic" zones with Ar « rA.9 ~ rsmOA<p would have been 
preferred, memory and computation time constraints re- 
quired the angular resolution be decreased by a factor of 
roughly two from the cubic value. The resolution of the 
three-dimensional runs was ultimately 360 x 360 x 1080 x 2 
{r X 8 X (f>, with the final factor of 2 representing the Yin 
and Yang sections of the grid), with rmax/f'min = 1.78. We 
deemed resolutions lower than this to be insufficient to cap- 
ture the structures expected in the runs with lower adiabatic 
indices. 

The instability seed for our simulations was a random 
multiplier between 0.95 and 1.05 applied independently to 
the density of each cell containing shocked ejecta; two- 
dimensional runs showed that the form of any initial per- 
turbation has a negligible effect on the final shape of the 
simulation remnant as long as the simulation is evolved long 
enough to include the instability saturation period discussed 
in the previous paragraph. 

Finally, ionization age of shocked ejecta was tracked to 
investigate its effect on observed morphology. As used in 
our simulations, the ionization time is a computationally 
inexpensive parameter that can then be used to provide a 
more accurate emission map than a basic density-squared 
model can. We first assumed that the electron density is 
proportional to the fiuid density throughout the remnant. 
During each cycle we calculated the temperature of every 
ejecta element; elements hotter than our cutoff temperature 
of 5 X 10^ K were assumed to have been shocked, and we 
updated the ionization time in each such cell by Tncw ~ 
Told + p ■ dt. Over the duration of our simulations, shocked 
ejecta never dropped below the cutoff temperature, so we 
could track shocked ejecta solely by temperature. We further 
assumed that iron is the sole contributor of free electrons to 
the ejecta; this is a very rough approximation, and future 
work might include additional species of ions for a more 
accurate picture. Using neon-like iron (16 free electrons per 
nucleus), the scale factor for ionization age is 

_ Pam ■ Na ■ 16 , 

"^scalc — t- 1~ 1^ ' 

55.8 

o It; , in9 2/3i-,-l/2 f Me \ ^^^ _g 

« 3.15x10 £5/ (^^^1 cm s. (7) 

where Na is Avogadro's number and 55.8 the molecular mass 
of iron. 

3.1 Dependence on time 

As mentioned in section 2, the only parameter affecting the 
morphology of a particular run is the scaled time. One of 
the primary questions addressed in this paper is whether the 
fleecy ejecta structures observed in type la SNRs are con- 
sistent with homogeneous ejecta. If the structures are solely 
due to hydrodynamic instabilities, their size and distribution 
could indicate the age of the remnant, providing a constraint 
on the three key parameters of the exponential model (Mcj , 
no, and -E51). For all three runs (7 — 5/3, 7 — 4/3, and 
7 = 6/5) we tracked the simulations longitudinally in time. 

Figure 1 shows the 7 = 5/3 run at three different times. 
Isosurfaces have been drawn of the three key interfaces - the 



Table 1. Interface Radii and Ratios from 3-D Data for 7 = 5/3 
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Here and throughout the paper, RcD is used to represent 
< RcT> >, the radius averaged over all Air steradians. 
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Figure 2. Interface locations as a function of time and adiabatic 
index for the exponential model in one dimension. The forward 
shock is tracked by thick solid lines, the reverse shock by thin solid 
lines, and the contact discontinuity by dashed lines. The curves 
for 7 = 5/3, 7 = 4/3, and 7 = 6/5 are in red, green, and blue 
respectively. The three curves for the contact discontinuity are 
separated by less than the width of the line used to show them. 

reverse shock, the contact discontinuity (where the ejecta 
fraction of fluid elements is 0.5), and the forward shock. 
All images in figure 1 are scaled to the same angular size, 
but the data presented in table 1 provide the relative scale. 
All images occur after the instabilities freeze out, so the 
production of new RT structures is negligible compared to 
the merging of existing structures. The three images have 
several distinguishing features, such as the proximity of the 
various fluid discontinuities and the structure at the forward 
and reverse shocks. 

The three interfaces consistently diverge from each 
other over the course of the simulation, illustrated by the 
one-dimensional simulations shown in figure 2: the contact 
discontinuity decelerates less than the reverse shock, and 
the forward shock decelerates less than either. Consequently, 
as time progresses there is less interaction between the two 
shock fronts and the Ray leigh- Taylor structures at the con- 
tact discontinuity. At the earliest time shown in figure 1, 
the interaction between the contact discontinuity and the re- 
verse shock is clear, as the features on the reverse shock are 
similar in angular size and spacing to the Rayleigh- Taylor 
structures just outside them. As the reverse shock recedes 
from the contact discontinuity, the shock front smooths out, 
the number of features and their radial extent shrinks, and 
their angular size grows. A similar process can be seen oc- 
curring at the forward shock: hy t — 2.0 the entire shock 
front has a radial dispersion of less than a single radial zone 
on the grid. 

The images in figure 1 show that the Rayleigh- Taylor 
structures in shocked ejecta trend toward larger angular 



© 2012 RAS, MNRAS 000, 1-17 



3-D morphology of la SNRs 5 




Figure 1. Isosurfaces for the reverse shock, contact discontinuity, and forward shock for the 7 = 5/3 run at three times: t = 0.12, 
t = 0.75, t = 2.0. The color scale for the contact discontinuity is absolute radial position, while the color scale for the forward and reverse 
shocks is percent difference from the average value for that interface as given in table 1. 
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Figure 3. Power spectra for the radial ejecta column density in 
the 7 = 5/3 simulation. The curves for early time {t ft: 0.15), 
middle time {t ^ 0.75) and late time (t = 2.0) are thin solid line, 
dashed line, and thick solid line respectively. The overlaid lines 
are best fits to power laws, P oc with all exponents ~ 3.9. 

All spectra have been normalized and smoothed. 
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Figure 4. Interface locations as a function of time and adiabatic 
index for the exponential model in one and three dimensions. 
Fluid discontinuities from 1-D runs are colored and marked as in 
figure 2. Interface locations in the 3-D runs are listed in tables 
elsewhere in the paper. The plots for 7 = 5/3, 7 = 4/3, and 
7 = 6/5 appear as the left, middle, and right panels in each triplet, 
respectively. Left: Close-ups of the 1-D interface locations around 
t = 0.15, with the 3-D interface locations plotted for comparison. 
The error bars show the maximum and minimum extent of the 
interface. To prevent overlap of 3-D error bars, some shocks are 
plotted slightly to the right of their proper position. Right: As 
before, but around t = 2.0; shifting here is to the left. 



size at later times. To further test this, we integrated the 
total shocked ejecta density in radial columns over all An 
steradians. We then used the SHTOOLS package^ to calcu- 
late power spectra for the generated image. The results are 
shown in figure 3, which plots the normalized power versus 
degree of spherical harmonic. The power spectrum's peak 
moves to lower values of I as the simulation progresses, from 
Z ~ 60 at the earliest time to Z ~ 30 at f = 2.0. The spectra 
between the peak and / ~ 200 can be fit to a power law in 
each case, with an exponent around 3.9. Higher than / ~ 200 
the bottleneck effect (Dobler et al. 2003) caused a steepen- 
ing of the power spectra in every case shown. Though the 
slope of the power spectrum does not substantially change 
with time, it is clear that power is transfered from higher 
wavenumbers to lower wavenumbers over the course of the 
simulation. 

All of the fits are significantly steeper than fitted power 
laws to the observed contact discontinuity of Tycho's SNR, 
P ~ k~^ '^ (Warren et al. 2005). Further, the difference is 
in the wrong direction: when the naturally two-dimensional 
surface of the CD is projected to a single line around the 
remnant, projection effects should act to smooth out some 
of the high-frequency power, steepening the power spectrum 
relative to that of the original 2-D surface. (See Appendix A 
for additional information.) While both power laws extend 
to a wavenumber of about 180, our peak occurs at a much 
higher wavenumber (/ ~ 30) than theirs (/ ~ 6). 

3.2 Dependence on adiabatic index 

Support abounds in the literature for the idea that efficient 
cosmic-ray acceleration impacts the evolution and morphol- 
ogy of SNRs. We find that increased shock compressibil- 
ity leads to dramatic differences in the morphology of the 

^ Available at http://shtools.ipgp.fr/ 
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shocked ejecta and of the two shock fronts. In this section we 
will describe the effect of the adiabatic index on four aspects 
of morphology: the shape and location of the contact discon- 
tinuity, the shape and location of the forward and reverse 
shock fronts, the power spectrum of the contact discontinu- 
ity, and the deceleration parameter. 

Under the assumption of spherical symmetry, reducing 
7 from 5/3 has the most immediate morphological effect of 
changing the shock locations of the remnant at any partic- 
ular time, as illustrated in figure 2. This does not apply to 
the contact discontinuity: the difference between the three 
1-D runs' contact discontinuities is less than the width of 
the line used to plot their position in figure 2. We explain 
this with an appeal to conservation of momentum. In the 
thin shell limit of Gull (1973), the entirety of the shocked 
ejecta and the swept-up ISM are contained at a single ra- 
dius. As the gas becomes less compressible, both the ejecta 
piston and the shocked ISM regions expand in volume rel- 
ative to the thin shell limit. This means that both shocked 
ejecta and the shocked ISM increase in mass, sweeping up 
additional matter that was beyond either shock front when 
the fiuid was more compressible. The additional inertia of 
the shell of ISM is roughly balanced out by the gained mo- 
mentum of the piston, and at a given time the radius of the 
contact discontinuity stays roughly constant as the compres- 
sion ratio drops from a = 00 to cr = 4. This effect is not 
permanent: the thickness of the regions of shocked fluid in- 
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Figure 5. Isosurfaces for the reverse shock, contact discontinuity, and forward slioclc at t = 2.0 for all three runs. The color scale for the 
contact discontinuity is absolute radial position, while the color scale for the forward and reverse shocks is percent difference from the 
average value for that interface as given in table 2. 
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creases constantly and reduces the appropriateness of the 
thin-shell approximation. After a few times t' (see section 
2), runs with lower compressibilities start to decelerate rela- 
tive to runs with greater compressibilities. The likely ages of 
Tycho's SNR and 1006 fall well short of the point of separa- 
tion, however, so for the purposes of this work the location of 
the contact discontinuity is independent of compressibility. 

Beyond spherical symmetry, the shape of the interfaces 
can also vary, in some cases substantially. Figure 5 shows 
all three three-dimensional runs at a scaled time of t = 2.0. 
As before, the remnants have been scaled to the same size, 
with the actual radii provided in table 2. 

In all of the tables in this paper, we list the location 
of the contact discontinuity as a single number, its average 
location. In actuality the contact discontinuity has a very 
complex shape that depends on both the remnant's age and 
the compressibility of the fluid, as can be seen in figures 1 
and 5. By t = 2.0, the CD is spread out over many radial 
zones (its radial extent is as high as 29% of its maximum ra- 
dius for the images in figure 5), and at many places around 
the remnant occurs multiple times in a single radial column 
(e.g., locations with RT mushroom caps above the base of 
the structure) . Using just the average location of the contact 
discontinuity is a great simplification. Reducing the struc- 
ture of the CD to a single number is nonetheless justified: as 
in one dimension (see figure 2), table 2 shows that the aver- 
age radius of the contact discontinuity in three dimensions 
is nearly constant across all three values of 7 at t = 2.0. 
Despite the independence of the CD's average radius and 
the compressibility of the gas, as 7 decreases there is a bias 
towards more ejecta close to the reverse shock; this is visible 
in figure 5 as a shift in color of the contact discontinuities 
away from yellow/white and toward red/black. 

The average location of the forward and reverse shocks 
in each of the 3-D runs is slightly less consistent with their 
radial positions in the corresponding 1-D simulations, as 
noted in figure 4 (the forward and reverse shocks have been 
plotted slightly to the left or right when needed for the sake 
of clarity). Both shock fronts are slightly more advanced in 
three dimensions - the reverse shocks are at a lower average 
radius, and the forward shocks at a greater average radius. 
The advanced location of the forward shocks is likely due to 
interaction between the shock front and Rayleigh- Taylor fin- 
gers at the contact discontinuity. As the adiabatic index de- 
creases and interaction between the two interfaces increases, 
protrusions appear at the forward shock that pull the aver- 
age location ahead of where it would be in one dimension. 
Figure 4 also shows error bars marking the maximum and 
minimum radii for each interface, further demonstrating the 
interaction and bubbles already mentioned. 

Figure 5 offers clear support for the effect of Rayleigh- 
Taylor fingers on the shape and location of the forward and 
reverse shocks. The 7 = 5/3 run shows very smooth for- 
ward and reverse shocks - the forward shock is located in 
the same radial zone everywhere in the remnant, while the 
reverse shock is spread over just a few zones. There is a 
large gap between both shocks and the contact disconti- 
nuity, the likely reason for the smoothness of both shock 
fronts. Although a major assumption behind our simula- 
tions was smooth ejecta, the second and third images in 
figure 5 demonstrate that smooth ejecta alone is not suffi- 
cient to guarantee smooth forward and reverse shocks. The 
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Figure 6. Normalized power spectra for the radial ejecta column 
density of the 3-D simulations at t = 2.0. The curves for 7 = 5/3, 
7 = 4/3, and 7 = 6/5 are in red, green, and blue respectively. 
The exponents to fitted power laws (not plotted here) are 3.9 for 
7 = 5/3, 4.0 for 7 = 4/3, and 3.5 for 7 = 6/5. 



7 = 4/3 run shows definite evidence of interaction at the 
forward shock: the majority of the interface is as smooth as 
its 7 = 5/3 counterpart, but a few tens of bumps can be seen 
where Rayleigh- Taylor fingers reached far enough outward 
to perturb the sphericity of the forward shock. The situation 
is repeated at the reverse shock, where features on the same 
angular scale as the RT structures are visible and the shock 
front itself has a radial extent of « 0.02 times its average 
radius in the data. When 7 = 6/5, the increased compress- 
ibility of the fluid has a marked effect on the shape of the 
remnant. There is abundant evidence of interaction between 
the contact discontinuity and the forward and reverse shocks 
(at the right edge of the image a Rayleigh- Taylor finger can 
even be seen in the process of creating one of the numerous 
bubbles visible on the forward shock). 

As the adiabatic index decreases from 5/3 to 1, the 
ejecta become more compressible and Rayleigh- Taylor struc- 
tures can be thinner. This results in more power at small 
wavelengths in the simulations with 7 < 5/3, illustrated 
in figure 6, which was created in the same manner as fig- 
ure 3. The figure shows an essentially monotonic increase 
in power at high wavenumbers as 7 decreases (the 7 = 4/3 
run's power spectrum peaks later, around I = 40, than do 
the power spectra of the other two runs, which peak around 
I = 30). Additionally, fitted power laws become shallower as 
the adiabatic index approaches 1, with the exponent drop- 
ping from 3.9 (7 = 5/3) to 3.5 (7 = 6/5). 

The deceleration parameter (m = vt/r) of the forward 
shock, as shown in figure 7, is another point of comparison 
between the radial structure in 1-D and its averaged equiv- 
alent in 3-D. Calculation of m in the 1-D case is simple, 
as one can easily track the motion of the forward shock, 
to grid-zone accuracy in space and arbitrary accuracy in 
time, and arrive at a value for the forward shock expansion 
velocity. For the 3-D runs, we used the Rankine-Hugoniot 
conditions for a strong shock to calculate the forward shock 
velocity in terms of the downstream pressure, upstream den- 
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Figure 7. The deceleration parameter as a function of time and 
adiabatic index for the exponential model. The curves for one 
dimension and 7 = 5/3, 7 = 4/3, or 7 = 6/5 are in red, green, 
or blue respectively. Overlaid are points showing the calculated 
deceleration parameter for the 3-D runs at selected times. 

sity, and adiabatic index. The various deceleration parame- 
ters are collected in figure 7. As in figure 4, the deceleration 
parameter in one dimension for each value of 7 is shown as a 
curve, with the corresponding 3-D values overlaid as points. 
At any particular time decreasing the adiabatic index re- 
duces the deceleration parameter, and a lowered adiabatic 
index reduces the time at which the remnant reaches a par- 
ticular value of the deceleration parameter. The curves and 
points of figure 7 are in close agreement everywhere, with 
the exception of the 7 = 6/5 run. There, the bubbles at the 
forward shock lead to a substantial solid angle where the 
shock normal isn't radial. This induced angle appears as a 
smaller radial expansion speed v, and so a smaller deceler- 
ation parameter. Restricting the averaging process to only 
those grid zones on the forward shock whose normal is within 
15° of the radial direction eliminates most of the bubbles but 
retains the largely spherical base visible in figure 5. It also 
moves the calculated 3-D deceleration parameters to within 
a few percent of the 1-D values, in line with the other values 
for 7; the corrected values are shown in figure 7, rather than 
the uncorrected numbers. 



4 OBSERVATIONAL IMPLICATIONS 

To facilitate comparison with the X-ray observations of Ty- 
cho and SN 1006, we now present line-of-sight projections 
from the three-dimensional data sets. We assume thermal 
emission from shocked ejecta to be proportional to the 
square of gas density. To mimic the effects of emission turn- 
on as shocked ejecta are ionized, a cutoff in ionization age is 
implemented by which emission from some ejecta elements 
can be excluded (see section 4.3). Synchrotron emission from 
the forward shock is visualized assuming that both relativis- 
tic electron energy density Ug and magnetic field energy den- 
sity ub are proportional to the pressure (i.e. nonlinear am- 
plification of the magnetic field). Then, since synchrotron 
volume emissivity is proportional to iteS'^"*'^'''^, with s 
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Figure 8. Images showing the effect that dynamical age has on 
the morphology of the remnant. The SE quadrant of SN 1006 
is included at the top left for comparison. All three remaining 
images come from the 7 = 5/3 run. Clockwise from the top right; 
t = 0.12, t = 0.75, t = 2.0. Image of SN 1006 taken from Cassam- 
Chenai et al. (2008). 

being the electron energy index (i.e. N{E) oc E~'') (Pachol- 
czyk 1970), we have cx pt^+^j/^. As s ^ 2.2 for both SN 
1006 and Tycho (Green 2009), the synchrotron emissivity is 
calculated as P^'*. X-ray synchrotron emission decays more 
rapidly away from the forward shock than radio emission, so 
our crude model results in a more diffuse shell of emission 
around the ejecta than would be present with a more refined 
treatment. In all images of projections in this section (with 
the exception of figure 10), ejecta emission is in white and 
synchrotron emission is in purple. 

Using these images, we take a more in-depth look at the 
distribution and character of emission from shocked ejecta, 
the shape and location of the projected contact discontinu- 
ity, and the strength of the limb brightening effect under 
ionization age cutoffs to emission. 

4.1 Fleece 

Figure 8 compares the SE quadrant of SN 1006 against pro- 
jected quadrants from the 7 = 5/3 run at the same three 
times shown (in 3-D) in figure 1. Immediately apparent is 
that the emission from shocked ejecta looks very similar to 
the fleecy complexes detected in X-rays in both Tycho and 
SN 1006. The earliest time shown in figure 8 (top right) is 
just after the end of the instability saturation period men- 
tioned in section 3. Ejecta structures at this time in the cen- 
tral region are visible as a filamentary network. In the time 
between the early image and the end of the simulation (bot- 
tom left of the figure) the Rayleigh- Taylor instabilities have 
grown, sheared, and merged into each other. Figure 9 com- 
pares the three different runs at t = 2.0 to the SW quadrant 
of Tycho's SNR to illustrate the effect of compressibility on 
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Figure 9. Images showing the effect that changing the adiabatic 
index 7 of the simulation has on the resultant remnant. The SW 
quadrant of Tycho's SNR is included at the top left for compari- 
son. Clockwise from the top right: the 7 = 5/3 run, the 7 = 4/3 
run, and the 7 = 6/5 run. All three projections are scaled to the 
correct relative size so interface locations can be directly com- 
pared. Image of Tycho taken from Warren et al. (2005). 

remnant morphology. The rounded mushroom caps of the 
Rayleigh- Taylor structures at 7 = 5/3 are still noticeable at 
7 = 4/3, though the center of the remnant is darker rela- 
tive to the limb brightening. By 7 = 6/6 the RT structures' 
caps are no longer generally round: the greater compressibil- 
ity has resulted in much longer, thinner fingers of shocked 
ejecta, and extensive interaction with the forward shock has 
bent or otherwise warped many of them. 

4.2 CD shape, relation to FS 

Figure 8 compared the 7 = 5/3 remnant at three different 
times to one quadrant of SN 1006. The structures visible 
at the center of SN 1006 are smaller in angular size than 
the RT structures at the center of the t = 2.0 image, but 
also less filamentary than the similarly-located ejecta at t = 
0.12. The RT fingers in both the 7 = 5/3 and the 7 = 
4/3 remnant are largely oriented along radial lines, as seem 
to be the structures at the edge of SN 1006. Furthermore, 
the structures at the edge of SN 1006 appear to be more 
discrete, that is, with a better-defined edge to each structure. 
The enhanced edges visible in the structures of SN 1006 
allow them to stand out against each other, as opposed to 
the less distinct haze outside the limb brightening in the 
t = 0.75 and t = 2.0 images. There is no obvious cutoff in 
limb brightening to mark a reverse shock in the image of SN 
1006, as the contrast is dominated by the large-scale gap in 
emission in the SE quadrant. 

As the compressibility of the fluid increases, interaction 
with the contact discontinuity can cause bubbles at the for- 
ward shock, noticeable in both the 7 = 4/3 and the 7 = 6/5 




Figure 10. Magnified views of two locations from the 7 = 6/5 
run where knots of ejecta seem to have overtaken the forward 
shock. 



images in figure 5. However, these bubbles are much fainter 
at their maximum radial extent than the shocked ISM at 
their base, so they do not generally show up in projection. 
One such bubble is visible in the 7 = 4/3 quadrant of fig- 
ure 9 as a thin bright rim of emission ahead of a darker patch. 
The situation is more extreme still with the 7 = 6/5 run, 
where the bubbles comprise a much larger fraction of the 
forward shock. Instead of the mostly smooth emission of the 
7 = 5/3 and 7 = 4/3 runs, the projected forward shock of 
the 7 = 6/5 run shows up as a chaotic network of projected 
bubbles, with single bubbles at the edge of the remnant too 
dim to appear in projection. It is for this run that the ap- 
proximation to X-ray synchrotron emission (discussed at the 
start of the section) is most telling; if the emission decayed 
more quickly, the diffuse shell of emission would resolve into 
a filamentary network tracing out the locations of the bub- 
bles at the forward shock. Comparison between figures 5 and 
9 implies that the locations around the rim of Tycho's SNR 
where emission from the forward shock is absent could be 
artifacts of interaction between the shocked ejecta and the 
forward shock. 

One feature captured by the 7 = 6/5 run, not present 
in the other two, is locations around the edge where knots 
of ejecta seem to be visible beyond the forward shock. Two 
examples have been enlarged and recolored for ease of view- 
ing in figure 10. The effect is an illusion caused by the faint 
bubbles of forward shock emission discussed in the previ- 
ous paragraph, and by projection effects; the contact dis- 
continuity is always inside the forward shock in the three- 
dimensional data. Even in the 7 = 6/5 run, such protrusions 
are rare, occurring only a few times around the rim of the 
remnant; they are absent entirely from the 7 = 4/3 and 
7 = 5/3 runs. With the smooth ejecta used in our runs, the 
ISM must be highly compressible before the forward shock 
is close enough to the contact discontinuity to be signifi- 
cantly affected by RT fingers. An alternative explanation 
was suggested by Orlando et al. (2012), who conclude that 
the separation between the CD and FS is an indication of 
ejecta structure rather than of cosmic ray acceleration: over- 
densities in the ejecta drive instability growth, resulting in 
more interaction at higher compressibilities than seen in our 
simulations. 

The knots pictured in figure 10 are roughly the same 
angular size, < 10°, as the protrusions around the rim of 
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Figure 11. Above, images of one simulated remnant with increas- 
ingly high cutoffs. Top left: an image with no cutoff in place; 
top right: r^in = 0.8; bottom right: Tmin = 1-6; bottom left: 
Tmin = 2.4. Using the scale factor from equation 7 and setting 
all parameters to unity, the three cutoffs correspond in physical 
units to Tmin = 2.52 X 10^ cm~^ s, Tmin = 5.04 X 10^ cm~^ s, and 
Tmin = 7.56 X 10^ cm~^ s respectively. All images are at t = 2.0 
for the 7 = 5/3 run. 

SN 1006 (in the SE, S, and SW) and Tycho (in the S and 
W). We conclude that these features can be generated by 
fluid instabihties alone, without any inhomogeneities present 
within the unshocked ejecta or the unshocked ISM. There 
are no features in the 7 = 6/5 run on the same scale as 
the shelf of thermal emission at the N rim of Tycho or the 
tiered structure in the NE polar cap of SN 1006, both of 
which are many tens of degrees across. The inability of the 
Rayleigh- Taylor instability to produce these features points 
at inhomogeneities in the ejecta, the magnetic field strength 
around the progenitor, or the ISM. 

4.3 Ionization age cutoffs, limb brightening, and 
angular variation 

Assuming that X-ray thermal emissivity is proportional to 
in visualizing data potentially overestimates emission 
from the remnant: not all matter is radiating in all wave- 
lengths at all times, due at least in part to deviations from 
ionization equilibrium. To account for this effect in the three- 
dimensional simulations, ionization age cutoffs were imple- 
mented during visualization, below which shocked ejecta 
were assumed to be X-ray faint. The effect on observed mor- 
phology is shown in figure 11, which (clockwise from the top 
left) sets the cutoff for emission at successively higher lev- 
els for the 7 = 5/3 run at t = 2.0. There is no cutoff for 
the top left quadrant (rmin = 0.0), and by the bottom left 
a high cutoff (rmin = 2.4) has eliminated about half of the 
shocked ejecta by volume. Since recently shocked ejecta is 
(for the 7 = 5/3 run) well within the innermost extent of 
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Figure 12. Illustration of the effect of a disrupted inner edge 
on the emission from a hollow shell. The solid red line shows a 
projection through a 2-D hollow shell; the radius of the inner 
edge is 0.75 times that of the outer edge. The dashed blue line is 
a projection through the same shell, but with the inner surface 
perturbed as 1 -|- 0.05 sin(4O0). The solid black line is similar, 
but with the perturbation given by 1 -|- 0.25 sin(406). Towards 
the center of the image the perturbations are nearly parallel to 
the line of projection, accounting for the noise visible in the two 
perturbed cases. 



the contact discontinuity, the cutoff preferentially eliminates 
smooth shells of recently shocked ejecta at low r and leads 
to raggedness at the inner edge of emission. 

This disruption of the inner boundary to emission pri- 
marily alters the strength of the limb brightening effect. This 
outcome is illustrated in figure 12, in which we generated 
two dimensional uniform hollow shells, perturbed the inte- 
rior edge, then projected them from two dimensions to one. 
Without any perturbations in place, the inner edge is clearly 
identifiable as the sharp peak in limb brightening; in pro- 
jections from three dimensions to two, a pronounced drop 
in brightness would occur. With a regular perturbation of 
amplitude 5% to the inner radius, the peak of emission no 
longer exactly traces the unperturbed reverse shock, even 
though the average radius is unchanged. Instead, the peak 
is outside that of the unperturbed case, and while the peak 
of emission is still plain the decay towards the center is no 
longer as pronounced. When perturbations are 25% of the 
average radius (an extreme case that isn't suggested by the 
three-dimensional data), the maximum of emission is now 
around that of the unperturbed case, but it is no longer 
a clear peak. Both the interior and exterior regions to the 
radius of maximum emission show gentle slopes, and the 
maximum itself is only ~ 150% of the emission from near 
the axis of symmetry of the shell. 

There is definite limb brightening occuring in Tycho's 
SNR, but the inner edge to the emission is indistinct. This 
effect is reproduced successfully by the lower two quadrants 
in figure 11, although the upper two quadrants (with no cut- 
off or a low one) still display the reverse shock as a localized 
drop in emission. Comparison to SN 1006 is less favorable; 
the limb brightening effect is either across too thin a re- 
gion (i.e. the northwest quadrant) or too broad a region 
(the southeast quadrant). 
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5 DISCUSSION 

This section elaborates on three topics that have been men- 
tioned only in passing so far. We examine the location of 
the interfaces (reverse shock, contact discontinuity, and for- 
ward shock) in the projected images of figures 8 and 9. We 
then combine available information throughout the paper 
with observed quantities to estimate the dynamical age of 
both SN 1006 and Tycho's SNR. Lastly, we discuss projec- 
tion effects and associated error with an eye toward future 
observations of remnants. 



5.1 Interface locations 

Addressing whether our remnants generate the radial struc- 
ture of Tycho's SNR and SN 1006 requires an analysis of our 
remnants mimicking the method of Warren et al. (2005) for 
locating interfaces. With only one "component" representing 
ejecta instead of nearly a dozen, however, we cannot repli- 
cate their procedure exactly. We first divided the remnant 
into 1440 angular wedges 0.25° in width, and created a ra- 
dial grid with 240 points along each angular wedge. Match- 
ing our remnant's radius to the 251" of Tycho's SNR, our 
radial resolution is roughly 1". This is finer than the value 
of 3" used in Warren et al. (2005) , but the resolution of our 
data is sufficient to allow it: Ar/r = 1.2 x 10""^, so the an- 
gular size of a single zone at the edge of the projection is 
just 0.3". 

With the projection partitioned into a grid, we located 
the contact discontinuity, forward shock, and reverse shock 
at each angular location. We identified the contact dis- 
continuity with the outermost radius where emission from 
shocked ejecta occurred. This method does not track the 
true contact discontinuity (shown in figures 1 and 5): as 
discussed in Appendix A, it biased toward protrusions at 
greater radius than the average value across the remnant. It 
is nonetheless very similar to the identification process used 
to find the CD in Warren et al. (2005) and Cassam-Chenai' 
et al. (2008). The forward shock was defined as the outer- 
most radius at which intensity reached half the maximum 
value for each radial spoke. This misses some of the filamen- 
tary structures that appear at lower values of 7, but does 
approximate the method of Warren et al. (2005). 

The reverse shock was found by treating the shocked 
ejecta as a hollow shell and identifying the maximum value 
along each radial line. This correctly identifies the location 
of the reverse shock for every combination of 7 and t except 
7 = 5/3, t = 2.0. In that instance, the density gradient in 
the shocked ejecta causes stronger emission near the contact 
discontinuity than near the reverse shock. This is visible in 
the upper left quadrant of figure 11 as a bright ring inside the 
edge of ejecta emission. The inner edge of emission becomes 
increasingly ragged or diffuse as 7 decreases or as we ex- 
clude freshly-shocked ejecta from emission calculations (see 
section 4.3 above). Any maximum value less than 50% of 
the overall maximum value was deemed to be unresolvable 
against the background of Rayleigh- Taylor fingers. Because 
of the Yin- Yang grid used for the simulations, overlap be- 
tween the two parts of the grid could generate an artificially 
high value for projected emission at the two angles where 
the overlap occurs. To prevent this artifact from affecting 
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Figure 13. Emission intensity as a function of radius for selected 
line-of-sight projections. All profiles have been scaled to the aver- 
age brightness of the central region (R < 0.2ijps) of that image. 
Top: intensity profiles for the 7 = 5/3 run. Peaking from left to 
right, t = 0.12 (dotted line), t = 0.75 (dashed line), and t = 2.0 
(solid line). The thin black lines trace out intensity due to shocked 
ISM. Bottom left: intensity profiles for the three runs a.t t = 2.0. 
Peaking from left to right, 7 = 5/3 (red line), 7 = 4/3 (green 
line), 7 = 6/5 (blue line). The thin lines trace out emissivity of 
shocked ISM. Bottom right: intensity profiles for the 7 = 5/3 run 
at t = 2.0, with ionization age cutoffs in place. Peaking from left 
to right, Tmin = 0.0 (thick solid line), Tmin = 0.8 (dotted line), 
Tmin = 1-6 (thin solid line), Tmin = 2.4 (dashed line). The thin 
black line is intensity of shocked ISM. 

the 50% threshhold of the reverse shock, the single highest 
value over all 1440 angular wedges was excluded. 

The angle-averaged brightness profiles for selected data 
sets are shown in figure 13, illustrating the effects of time, 
compressibility, and an ionization age cutoff on emission. 
The effects of the ejecta gradient mentioned previously are 
visible as the plateau in the late-time curve in the top panel. 

The radii for the three fluid discontinuities, averaged 
around the rim of each projection, are gathered in table 3, 
based off of the preceding analysis of the line-of-sight projec- 
tions. In this table and for the rest of the paper, we use the 
notation of Cassam-Chenai et al. (2008), in which -Rcd de- 
notes the projected radius of the contact discontinuity and 
-Rcd the average location of the discontinuity in three di- 
mensions. Comparing Rys between the runs with 7 = 5/3 
and 7 = 4/3, it is clear that increasing the compressibility 
of the fluid has a significant effect on the overall size of the 
remnant at similar dynamical times, just as with the one- 
dimensional runs presented in figure 2. The forward shock 
decelerates less than the reverse shock and the contact dis- 
continuity as the remnant evolves, and all ratios in Table 3 
drop with time. Lowering the adiabatic index has the oppo- 
site effect, bringing the interfaces closer together and raising 
the ratios in table 3. 

From shock positions alone it may be difficult to distin- 
guish between a younger remnant with less efficient cosmic 
ray acceleration and a more evolved remnant with more ef- 
ficient particle acceleration. At roughly the age of Tycho's 
remnant, the ratio Rco/Rfs for the 7 = 5/3 run is lower 
than observed in Tycho, as is Rns/Rps- With the 7 = 4/3 
simulation, Rcd/Rfs is close to that observed in Tycho, 
but the reverse shock is too close to the forward shock. If 
the forward shock were more efficient at cosmic ray accel- 
eration than the reverse shock, the effective adiabatic index 
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Table 3. Interface radii and ratios, from 2-D projections and 3-D data 
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0.892 


<1% 
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1.391 


1.616 


0.624 


0.861 


1.005 


1.229 
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0.351 
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1% 
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0.946 
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6% 
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0.877 
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9% 
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1.393 


1.449 


0.769 


0.961 


1.112 


1.224 


1.438 


0.773 


0.851 
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of shocked ejecta and shocked ISM would differ, leading to 
varying compressibilities ahead of and behind the contact 
discontinuity. Previous work offers support for such a sit- 
uation. Decourchelle, Ellison & Ballet (2000) provided evi- 
dence that low magnetic fields lead to negligible acceleration 
of cosmic rays. In the absence of magnetic field amplifica- 
tion by the reverse shock, the frozen-in magnetic field of the 
ejecta should be attentuated by expansion, and therefore 7efi 
should be approximately 5/3. Ellison, Decourchelle & Ballet 
(2005) elaborated on the bijective relationship between mag- 
netic field amplification and DSA, and while studies have 
been done of magnetic field amplification at Tycho's for- 
ward shock (Volk, Berezhko & Ksenofontov 2005), we are 
unaware of observations suggesting magnetic field amplifi- 
cation at Tycho's reverse shock; Hughes, Rakowski & De- 
courchelle (2000) did, however, find evidence of efficient cos- 
mic ray acceleration at the reverse shock of IE 0102.2-7219. 
Our conclusion echoes that of Warren et al. (2005), which 
found little evidence for particle acceleration at the reverse 
shock of Tycho by comparing 1-D hydrodynamic models to 
the discontinuity radii deterimined in that paper. 



5.2 Implications for Tycho and SN 1006 

The evolution of the exponential model is governed entirely 
by the scaled age t, so attempts to compare the model to ob- 
servations must determine the dynamical age of the remnant 
in question. In the previous sections we have examined the 
effect of scaled time and adiabatic index on the appearance 
and structure of our simulated remnants, including interface 
locations, power spectra, and deceleration parameters. Now 
we will combine these results with observed quantities to 
estimate the dynamical ages of Tycho and SN 1006, and to 
discuss the efficiency of cosmic ray acceleration in both rem- 
nants. We consider the fleecy ejecta structures in the SNRs, 
dynamical quantities of the ejecta, and the average radial 
locations of the fluid discontinuities in determining the age 
of the two SNRs. We also use deceleration parameters in or- 
der to constrain the effective adiabatic index at the forward 
shock. 

Tycho's SNR was first observed in 1572, giving it (at 
time of writing) a real age of 439 years. This corresponds to a 
scaled age of t = 439 yr/T' « 1.77 {Me/Mch)~'^^^Eli^nl^'-^. 
SN 1006, with a real age of 1005 years, has a scaled age 
of t = 4.05 {Me/Mch)~^^'^ El(^nl^'\ Despite the difference 
in age, our simulations strongly suggest that Tycho be dy- 
namically older than SN 1006. The fleecy structures present 



in Tycho are, on average, larger in angular size than those 
observed in SN 1006. In our sinmlated remnants and their 
projections, the angular size of the fleecy ejecta structures is 
almost completely determined by the dynamical age of the 
remnant (see flgures 1, 3, and 8); the size of the structures 
depends only weakly, if at all, on the compressibility of the 
fluid (figures 5 and 6). 

One way to approach the dynamical age of the remnant 
is consideration of the motion of the ejecta. Three quantities 
are relevant here: (i) the free expansion velocity of unshocked 
ejecta just ahead of the reverse shock, (ii) the bulk radial 
motion of shocked ejecta near the reverse shock, and (ill) 
the bulk motion of the ejecta at the contact discontinuity. 

(i) The expansion velocity at the reverse shock of SN 1006 
is 7026 ± 13 km s'^iv = OSSiM^/MchY^^E^^^^) based 
on Doppler-shiftcd absorption lines through the remnant 
(Hamilton et al. 2007) (the shape of the lines points to asym- 
metry in the remnant, however, so it is uncertain if this ve- 
locity is representative of the reverse shock everywhere). The 
shock is unlikely to be accelerating cosmic rays efficiently, 
for reasons outlined in the previous subsection; the effective 
adiabatic index should therefore be close to 7 = 5/3. In 
this case, the location of the reverse shock is consistent in 
one dimension and in three. From 1-D runs, then, we find 
that a free expansion velocity of 0.83 at the reverse shock 
corresponds to a dynamical age of f = 1.0 for SN 1006 for 
canonical values oi E^i = 1 and Me = Mch- 

(ii) No observations like those of Hamilton et al. (2007) 
have been reported for Tycho, but Hayato et al. (2010) fitted 
pairs of red- and blue-shifted absorption lines to the Fe Ka 
spectral feature of Tycho, and measured the radial velocity 
by quantifying the relative shift between the center of the 
remnant and the edge. They report expansion velocities of 
4000 ±300 km s'^v = QA7 [Me/ Mchf^'^ E^-^''^) for Fe Ka, 
which should be associated with freshly shocked ejecta at 
the reverse shock. The velocity of recently shocked ejecta in 
our 7 = 5/3 simulation at t = 0.75 {t = 2.0) is v = 0.58 {v = 
0.21). Interpolation between the two suggests a dynamical 
age of 1.1 for Tycho. 

(iii) Using a similar process for the Si He/3 feature, 
Hayato et al. (2010) found the expansion velocity to be 
4700 ± 100 km s"i(t) = Q.56{M^/Mci,Y^'^ E~l^^) for sili- 
con, which they associate with the contact discontinuity. In 
1-D the radial velocity of the contact discontinuity is 0.56 at 
t = 0.72. However, in multiple dimensions fluid instabilities 
cause a range of velocities and densities in the mixing region. 
The greatest emissivities should occur in the forward tips 
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Table 4. SNR 1006 parameters, assuming Me/Mch = 1 



Me/Mch 




(1) 


t'- > 


no (cm ■^J'--^-' 




1.0 


1.0 


0.83 


0.98 


0.014 


8.84 


1.0 


1.5 


0.68 


1.35 


0.020 


9.31 


l.U 


o.u 


0.48 


2.15 




10 42 


1.0 


4.0 


0.42 


2.50 


0.030 


11.03 


1.0 


6.0 


0.34 


3.18 


0.031 


12.04 


1.0 


9.0 


0.28 


3.67 


0.028 


13.37 


1.5 


2.0 


0.72 


1.25 


0.029 


9.13 


1.5 


2.5 


0.64 


1.48 


0.035 


9.39 


1.5 


3.0 


0.59 


1.66 


0.037 


9.69 


1.5 


4.0 


0.51 


2.00 


0.042 


10.17 


1.5 


6.0 


0.42 


2.50 


0.043 


11.03 


1.5 


9.0 


0.34 


3.10 


0.047 


11.99 


(1) 7026 km 


in 


scaled units. 







(2) Time at which reverse shock velocity equals vns- 

(3) Density required for the scaled time t to correspond to 1001 

yr. 

(4) At scaled time t. 

of Rayleigh- Taylor structures, which expand more rapidly 
than the rest of the contact discontinuity. The highest ion- 
ization ages should also occur in the R-T structures. In the 
7 = 5/3 run, the 10 per cent of the ejecta with the highest 
ionization age has a radial velocity of 0.65 aX t = 0.75; this 
decreases to v — 0.26 hy t = 2.0. Interpolating between the 
two numbers, the results of Hayato et al. (2010) applied at 
the contact discontinuity of our 3-D simulations imply that 
Tycho's dynamical age is 1.0. 

The estimated age around 1.0 for both Tycho and SN 1006 
conflicts with the size of the ejecta structures, which should 
be a marker for relative age. The reason for this difference 
is not clear. 

We now assess the ramifications of observed fluid ve- 
locities on the parameters (no. Me, and £51) of the two 
SNRs. As mentioned in section 2, the scaling factor for time 
depends inversely on both the explosion energy and the in- 
terstellar density (it also depends on the ejecta mass, but 
we see little reason why this should deviate from Mch)- 
There is a good deal of evidence that the ISM around Tycho 
is denser than that around SN 1006: the inferred densities 
around Tycho from X-ray (Cassani-Chcnai' ct al. 2007), 7- 
ray (Volk, Berezhko & Ksenofontov 2008), and optical ob- 
servations (Kirshner, Winkler & Chevalier 1987) all support 
a value for no of 0.1-0.3 cm~^ in the west, with higher densi- 
ties in the east. Assuming £51 = 1, the scaled and physical 
ages of Tycho match if no = 0.18 cm~^, consistent with 
observations. 

Around SN 1006, observations suggest a lower interstel- 
lar density. Acero, Ballet & Docourcholle (2007) found the 
X-ray emission measure of the forward shock and calculated 
an ISM density of 0.05 cm~^ in the southeast quadrant, 
and argued that the density was roughly constant every- 
where except the filamentary northwest rim. Using expan- 
sion data, Katsuda ct al. (2009) determined the density at 
the northeast to be 0.085 cm~^. A very low ISM density, 
no = 0.014 cm^^ if E51 — 1, is required to match the dy- 
namical and physical ages (see table 4). This is lower than 
either observational estimate, though it is within model- 
dependent uncertainties of the value presented by Acero, 



Ballet & Decourchelle (2007). Also visible in table 4 is that 
the physical age of SN 1006 cannot be matched to the veloc- 
ity measurements of Hamilton et al. (2007) for any ambient 
density greater than about 0.030 cm~^, which itself requires 
t > 2.15 and E51 > 3.0. We can impose an additional con- 
straint on the parameters of SN 1006 by considering its size: 
at a distance of 2.18 ± 0.08 kpc (Winkler, Gupta & Long 
2003), and with an average radius of 14.5' in the SE quadrant 
(Cassam-Chenai et al. 2008), the radius of the forward shock 
is 9.19 parsecs. Under the assumption that Me/Mch = 1, the 
explosion energy £^51 = 1.4 and ISM density no — 0.019 are 
consistent with SN 1006's age, reverse shock velocity, and 
size. These parameters suggest a scaled age of 1.3 for SN 
1006. If we relax the assumption that Me = Mch, we again 
find a consistent solution. Table 4 shows the same calcula- 
tions if Me/Mch = I.5. It can bo seen that the upper limit 
on no is 0.047 cm^"^ in this case. The parameters that match 
the age, size, and reverse shock velocity are E^i =2.1 and 
no = 0.030 cm~^; they also imply that SN 1006's dynamical 
age is 1.3. Such a massive and energetic supernova could be 
consistent with a white dwarf merger origin for SN 1006. 

The deceleration parameter, like the separation between 
fluid discontinuities and the free expansion velocity of ejecta, 
is a monotonic function of time. As evidenced by figure 7, 
it is also sensitive to the effective adiabatic index: for any 
particular time, decreasing the adiabatic index also lowers 
the deceleration parameter. Tycho's deceleration parameter 
is higher along the NW-SW rim (where the remnant is close 
to spherical symmetry), 0.59 ± 0.12 (Katsuda et al. 2010), 
than is SN 1006's at the NE rim (the only quadrant where 
such a study has been performed in X-rays), 0.54 ± 0.06 
(Katsuda et al. 2009), although the two numbers are sepa- 
rated by a single standard deviation. Since the deceleration 
parameter should decay from 1.0 to 0.4 as a remnant transi- 
tions from free expansion to the Sedov phase, these numbers 
suggest that Tycho is slightly dynamically younger than SN 
1006. The large uncertainties in both values do leave open 
the possibility that the situation is reversed - that the cor- 
rect deceleration parameter for Tycho is lower than that for 
SN 1006 - but the data presented in figure 7 offer another 
interpretation consistent with observations. We posit that, 
if SN 1006 is dynamically younger than Tycho - or even 
the same dynamical age - it must be more efficiently accel- 
erating cosmic rays at its forward shock, lowering both its 
deceleration parameter and effective adiabatic index below 
Tycho's. 

To quantify the difference in adiabatic index of the two 
remnants, we turn to the forward shock and its separation 
from the contact discontinuity. There is substantial evidence 
that both remnants are efficiently accelerating cosmic rays 
at their forward shock, so we exclude the 7 = 5/3 run and 
focus on the 7 = 4/3 and the 7 = 6/5 runs. At i = 1.0, 
RcD ■ Rfs is 0.96 and 0.99, respectively, for 7 = 4/3 and 
7 = 6/5. This ratio for SN 1006 is 0.98 in the NE polar cap 
(Cassam-Chenai et al. 2008). The effective adiabatic index 
in SN 1006 should therefore be just over 1.2 - at the very 
least, nmch closer to 7 = 6/5 than to 7 = 4/3. The esti- 
mated adiabatic index is less than 6/5 if the age estimate 
is increased to 1.3, as suggested by Table 4. In the case of 
Tycho, RcD : Rfs along the western rim (where the rem- 
nant appears mostly spherical) is around 0.95. We expect 
an effective adiabatic index of just over 4/3 for Tycho (this 
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number increases if the dynamical age decreases, and vice 
versa). We note that these numbers agree broadly with the 
results of Kosenko, Blinnikov & Vink (2011), who used 1-D 
simulations to calculate a range of allowable compressibili- 
ties that could match the radial morphology of Tycho and 
SN 1006, and found that the upper end of SN 1006's range 
was higher than the upper end of Tycho's. 

As mentioned in section 3.1 and tabic 3, the reverse 
shock and contact discontinuity are constantly receding rel- 
ative to the forward shock; the separation between the in- 
terfaces could therefore serve as a probe of dynamical age. 
No reverse shock has been directly observed in SN 1006, 
but Warren et al. (2005) identified Fe Ka emission with 
the reverse shock in Tycho. Since the location of the con- 
tact discontinuity is far less dependent on the adiabatic 
index than are the locations of the forward and reverse 
shocks, _Rr,s : -Rcd requires fewer assumptions as an in- 
dicator of Tycho's dynamical age than docs the more com- 
monly quoted figure of .Rrs : .Rfs- For Tycho, this value is 
0.73/0.96 = 0.76 (Warren et al. 2005). The same ratio for 
the 7 = 5/3 run a.tt = 0.75 is 0.83, and has dropped to 0.73 
hy t = 2.0. Interpolation between the two yields a dynami- 
cal age for Tycho of 1.6. The ratio /?rs : Rps can be used 
to estimate Tycho's dynamical age, but requires an assump- 
tion about the compressibility of the ejecta and ISM. Using 
7 = 5/3 for the ejecta (and Rks) and 7 = 4/3 for the ISM 
(and Rys), the appropriate ratio is 0.82 at t = 0.75 and 0.68 
&tt = 2.0. The ratio in Tycho is 0.73 (Warren et al. 2005), 
which occurs at t = 1.6. This result does not change signif- 
icantly if more compressible ISM is assumed. Intriguingly, 
velocity measurements and radial morphology each give a 
consistent estimate for the age of Tycho's remnant, but the 
two estimates do not agree with each other. 

Taken together, the results discussed in this section 
paint the following scenario for the relative dynamical ages 
and effective adiabatic indices of the remnants of Tycho's 
SN and SN 1006. The size of ejecta stuctures and Tycho's 
radial morphology suggests that it is the older remnant, but 
available velocity information implies that both remnants 
have roughly equal dynamical ages. In either case, SN 1006 
is more efficiently accelerating cosmic rays at its forward 
shock. We calculate a dynamical age for SN 1006 of 1.3 based 
on the free expansion velocity of unshocked ejecta and on 
the known size of the remnant. At this age, the separation 
between the forward shock and the contact discontinuity im- 
plies an effective adiabatic index of 6/5. If the scaled age of 
Tycho's SNR is 1.0 as suggested by expansion velocity of its 
ejecta, we find that 7efi « 4/3 at its forward shock. Given 
the simplicity of our model, the disparate age estimates in 
the case of Tycho, and the numerous factors affecting the 
remnants' actual expansions, however, further study is war- 
ranted before firm conclusions can be drawn. 



6 CONCLUSIONS 

We have performed high resolution three-dimensional sim- 
ulations of a Type la supernova remnant using an expo- 
nential ejecta profile and assuming homogenous ejecta and 
ISM. The Rayleigh- Taylor instability is clearly capable of 
generating the fieecy structures observed in Tycho's SNR. 
Comparison against SN 1006 is qualitatively less favorable 



due to the orientation and spacing of the structures in that 
remnant. The ejecta and ISM in all simulations were smooth 
at initialization, implying that dumpiness is not a necessary 
condition to generate the structures observed in Tycho. This 
result depends critically on evolving the simulations long 
enough time for the instabilities to saturate, independently 
of the compressibility of the fluid. While our simulations re- 
produced the central regions of both remnants well, they 
require a relatively high cutoff in ionization age to capture 
the indistinct limb brightening seen in Tycho, and qualita- 
tively fail to match SN 1006. 

After consideration of several observables tied to dy- 
namical age - such as expansion rate of ejecta at the reverse 
shock, average radii of fluid discontinuities, and deceleration 
parameter for select regions around the rim - we find that 
observed parameters are inconsistent in the case of Tycho's 
SNR, with radial structure and ejecta morphology present- 
ing a different age (i = 1.6) than observed ejecta velocities 
{t = 1.0). Taking the lower value as the more accurate, Ty- 
cho's effective adiabatic index is slightly higher than 4/3 
at its forward shock. SN 1006's dynamical age is approxi- 
mately 1.3, and accurate measurements of its distance and 
size allow us to calculate E51 = 1.4 and no = 0.019 cm~"^. 
Acceleration of cosmic rays at the forward shock of SN 1006 
is more efficient than at that of Tycho's SNR: 7cff is 6/5. 

Although shape of the contact discontiimity depends on 
adiabatic index and compressibility of the ejecta, its aver- 
age radius does not for the period of time covered by our 
simulations. Further, the shape of the CD, especially given 
the close proximity to the forward and reverse shocks at low 
7, affects the shape of the shock fronts. For the lowest value 
of 7 used, knots of ejecta appear to protrude outside the 
forward shock, but in reality lie just inside faint bubbles of 
emission from shocked ISM. This does not happen on an- 
gular scales as large as the shelf of emission in NE Tycho 
or the polar region of emission in NE SN 1006. Since insta- 
bilities and smooth ejecta cannot generate those features, 
something else (inhomogeneity in ISM, ambient magnetic 
fields, or asymmetric explosion) must be necessary. Studies 
including inhomogenous ejecta or an ambient magnetic field 
to azimuthally affect cosmic ray production offer additional 
insight into the laxge-scale structure of both of these rem- 
nants. 
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APPENDIX A: DIFFERENCES BETWEEN 2-D 
PROJECTIONS AND 3-D DATA 

Projections of an irregular three-dimensional spherical struc- 
ture are biased towards protrusions from the average ra- 
dius, so the average radius of the contact discontinuity 
differs from the observed radius in projection. This bias 
can be parameterized by a correction factor ^proj, where 
RpTO] = (1 + Cproj)J?true. The corrcctiou factor depends 
strongly on the degree of structure present in a surface, 
so the forward and reverse shocks of our model remnants 
should see minimal correction compared to the contact dis- 
continuity. In Warren et al. (2005) several factors (e.g. the 
observed power spectrum and expected length of Rayleigh- 
Taylor fingers) were considered in determining a correction 
factor of 6% for Tycho, and resulted in a revision from ob- 
served ratios (1 : 0.96 : 0.73) to "true" ratios (1 : 0.93 : 0.71). 
Drawing on the work of Dwarkadas & Chevalier (1998) and 
Wang & Chevaher (2001), Cassam-Chenai et al. (2008) ar- 
rived at a projection correcting factor of 10% for SN 1006. 
With access to both projections and fully three-dimensional 
data, we can directly evaluate the errors of this method. 

The projection correcting factor is related to the pro- 
jected and true discontinuity locations by {Rcyi/Rfs) = 
(RcD / Rfs) • (1 + C); ^ similar equation exists for the re- 
verse shock. Table 3, in addition to the interface locations 
in the projections, also lists both the true discontinuity loca- 
tions and the projection correction factor for each 'y/t pair. 
Both the CD/FS and RS/FS ratios arc higher in projection 
than they are in 3-D for every simulation at almost every 
time; the two exceptions for 7 = 6/5 are likely caused by 
the filamentary nature of emission at the forward shock and 
the extensive interaction between the fluid discontinuities. 
Additionally, ^ is very nearly for the reverse shocks. This 
result is in keeping with the relative smoothness of both the 
forward and reverse shocks as seen in figures 1 and 5. War- 
ren et al. (2005) predicted a ^proj of 3% for the reverse shock 
of Tycho, an overestimation of the structure present in the 
reverse shock. 

For the contact discontinuity, the correction factors 
range from 6% at the earliest times to 12-13% at the latest, 
showing remarkable consistency across 7 despite the changes 
to the shape of both the forward shock and the contact dis- 
continuity. Cassam-Chenai et al. (2008) suggested that a ^ 
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upwards of 10% would be necessary to match observations of 
SN 1006 (7?cD /Rfs = 0.96 for the SE quadrant where mini- 
mal nonthermal emission is seen) with their onc-dimensional 
simulations. Their runs were carried out to a scaled time of 
t w 1.2, depending inversely on the interstellar density in 
the environment of SN 1006. Given the lack of synchrotron 
emission observed over this region, for our hydrodynamic 
models to match that shock ratio we would require an ex- 
ceptionally low ambient density in those regions (reducing 
the dynamic age of the remnant) in conjunction with effi- 
cient acceleration. 
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